Libraries to load:

library(dplyr)
library(INLA)
library(inlabru) 
library(sf)
library(terra)


# load some libraries to generate nice map plots
library(scico)
library(ggplot2)
library(patchwork)
library(mapview)
library(tidyterra)

The data

In this practical, we will revisit the data on the Pacific Cod (Gadus macrocephalus) from a trawl survey in Queen Charlotte Sound. The pcod dataset is available from the sdmTMB package and contains the presence/absence records of the Pacific Cod during each surveys along with the biomass density of Pacific cod in the area swept (kg/Km\(^2\)). The qcs_grid data contain the depth values stored as \(2\times 2\) km grid for Queen Charlotte Sound.

The dataset contains presence/absence data from 2003 to 2017. Lets consider year 2003 for now.

library(sdmTMB)

pcod_df = sdmTMB::pcod 
qcs_grid = sdmTMB::qcs_grid

Then, we create a sf object and assign the right coordinate reference to it:

pcod_sf =   st_as_sf(pcod_df, coords = c("lon","lat"), crs = 4326)
pcod_sf = st_transform(pcod_sf,
          crs = "+proj=utm +zone=9 +datum=WGS84 +no_defs +type=crs +units=km" )

We convert the covariates into a raster and assign the same coordinate reference:

depth_r <- rast(qcs_grid, type = "xyz")
crs(depth_r) <- crs(pcod_sf)

Spatio-temporal LGOCV comparison

Model fitting

Now lets compare two different space-time models using LGOCV and some information criteria metrics. The general model structure is given by:

\[ \begin{aligned} y(s,t)|\eta(s,t)&\sim\text{Binom}(1, p(s,t))\\ \eta(s,t) &= \text{logit}(p(s,t)) \\ \end{aligned} \]

  1. Model 1 - time iid effect We consider a separable space-time model with a linear predictors given by:

\[ \eta(s,t) = \beta_0 + f_1(\text{depth}(s)) + f_2(t) + \omega(s) \]

  • \(f_1(\text{depth}(s))\) is a smooth covariate effect of depth (modeled using a RW2 model)

  • \(f_2(t)\) is an IID effect of time

  • \(\omega(s)\) is Matérn random field.

The first step is to define the mesh and the spde model

Construct the mesh and the SPDE model

mesh = fm_mesh_2d(loc = pcod_sf,    
                  cutoff = 1,
                  max.edge = c(10,20),     
                  offset = c(5,50),
                  crs = st_crs(pcod_df))   

spde_model =  inla.spde2.pcmatern(mesh,
                                  prior.sigma = c(1, 0.5),
                                  prior.range = c(100, 0.5))

create time index and the grouped variable

To use the RW2 model the covariate has to be groupes:

depth_r$depth_group = inla.group(values(depth_r$depth_scaled))

we also define a time index from 1 to in the data frame

pcod_sf = pcod_sf %>%
     mutate(time_idx = match(year, c(2003, 2004, 2005, 2007, 2009, 2011, 2013, 2015, 2017)),
         id = 1:nrow(.)) # Observation id for CV

##Implement Model 1##

Warning{{< bi pencil-square color=#c8793c >}} Task

Implement the model in inlabru

  1. Define the components
  2. Define the formula
  3. Define the likelihood model using the bru_obs() function
  4. Run the model
# Model components
cmp_spat = ~ Intercept(1) + 
  covariate(depth_r$depth_group, model = "rw2", scale.model = TRUE)+
  trend(time_idx, model = "iid")+
  space(geometry, model = spde_model)

# Linear predictor
formula_spat = present ~ Intercept  + trend + space + covariate

# Observational model
lik_spat = bru_obs(formula = formula_spat, 
              data = pcod_sf, 
              family = "binomial")

# Fit Model 
fit_spat = bru(cmp_spat,lik_spat)
Note

Note that there are some survey locations in certain years that fall outside the depth raster region. inlabru will input these missing covariate values using the nearest available value. This can be computationally expensive, but you can avoid it by supplying a raster layer that encompasses all of your data points (e.g., by pre-imputing these missing values with your preferred method of choice).

One way of doing this in the inlabru framework is to use the bru_fill_missing() function:

# Select the raster of interest
depth_orig = depth_r$depth_group
re <- extend(depth_orig, ext(depth_orig)*1.05)
# Convert to an sf spatial object
re_df <- re %>% stars::st_as_stars() %>%  st_as_sf(na.rm=F)
# fill in missing values using the original raster 
re_df$depth_group =  bru_fill_missing(depth_orig,re_df,re_df$depth_group)
# rasterize
depth_filled <- stars::st_rasterize(re_df) %>% rast()
plot(depth_filled)

We have now fit the model and want to check the results.

Warning{{< bi pencil-square color=#c8793c >}} Task

Is the effect of depth significant? Does it seem important to have a non-linear model?

Inspect the estimated time effect and

Use the predict function to inspect the estimated spatial effect.

Use the depth_r raster to define the points in space where to predict

pxl = st_as_sf(data.frame(crds(depth_r)), coords = c("x","y") ,
               crs  = st_crs(pcod_sf))
# covariate effect

p_cov = fit_spat$summary.random$covariate %>%
  ggplot() + geom_ribbon(aes(ID, ymin = `0.025quant`, ymax = `0.975quant` )) +
  geom_line(aes(ID,mean)) + ggtitle("Covariate effect")

# time effect
p_time = fit_spat$summary.random$trend %>%
  ggplot() + geom_errorbar(aes(ID, ymin = `0.025quant`, ymax = `0.975quant` )) +
  geom_point(aes(ID,mean)) + ggtitle("Time effect")

#space effect
pred_space = predict(fit_spat, pxl, ~ space)

p_space_mean = ggplot() + gg(pred_space, aes(color = mean))
p_space_sd = ggplot() + gg(pred_space, aes(color = sd))
  1. Model 2 - spatiotemporal field We consider a separable space time model with a linear predictor given by:

\[ \eta(s,t) = \beta_0 + f_1(\text{depth}(s)) + \omega(s,t) \]

  • \(f_1(\text{depth}(s))\) is a smooth covariate effect of depth (RW2)
  • \(\omega(s,t)\) is a space-time Matérn spatial field with AR1 time component \[ \omega(s,t) = \phi\ \omega(s,t-1) + \epsilon(s),\qquad \epsilon(s)\sim\text{GF}(\sigma_{\epsilon},\rho_{\epsilon}) \]
Warning{{< bi pencil-square color=#c8793c >}} Task

Implement this second model in inlabru

  1. Define the components For the AR1 model use this following PC prior for the correlation parameter \(\phi\)
# PC prior for AR(1) correlation parameter
h.spec <- list(rho = list(prior = 'pc.cor0', param = c(0.5, 0.1)))
  1. Define the formula
  2. Define the likelihood model using the bru_obs() function
  3. Run the model (This model can take a couple of minutes to run)
# Model components
cmp_spat_ar1 = ~ Intercept(1) + 
  covariate(depth_filled$depth_group, model = "rw2", scale.model = TRUE)+
  space_time(geometry,
        group = time_idx ,
        model = spde_model,
        control.group = list(model = 'ar1',hyper = h.spec))

# Linear predictor
formula_spat_ar1 = present ~ Intercept  + covariate  + space_time

# Observational model
lik_spat_ar1 = bru_obs(formula = formula_spat_ar1, 
              data = pcod_sf, 
              family = "binomial")

# Fit Model 
fit_spat_ar1 = bru(cmp_spat_ar1,lik_spat_ar1)